ACCESS overestimates concentration in polar region while underestimates in
import cosima_cookbook as cc
from cosima_cookbook import explore
import numpy as np
import dask.array as da
import xarray as xr
import pandas as pd
import os.path
import glob
from tqdm.notebook import tqdm
#------------------------
#----------------------
import logging
logging.captureWarnings(True)
logging.getLogger('py.warnings').setLevel(logging.ERROR)
from dask.distributed import Client
# change to your own directory on /g/data/:
figdir = '/g/data/hh5/tmp/rsp599_tmp/ACESSOM201/figures/BGC_IAF/'
Recommended plotting libraries
%matplotlib inline
%config InlineBackend.figure_format='retina'
#Packages for plotting
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.dates as mdates
import matplotlib.colors as mcolors
from matplotlib.ticker import AutoMinorLocator
import matplotlib.cm as mcm
import matplotlib.gridspec as gridspec
#
import IPython.display
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import geopandas
#
land_50m = cft.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='black', facecolor='gray', linewidth=0.5)
client = Client()
client
# client = Client(processes=2, threads_per_worker=1,
# n_workers=2, memory_limit='10GB')
# client
Client-6ec8b31e-86a2-11ed-a0db-fa163eefa180
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: /proxy/42597/status |
eb9aa26a
| Dashboard: /proxy/42597/status | Workers: 4 |
| Total threads: 16 | Total memory: 44.92 GiB |
| Status: running | Using processes: True |
Scheduler-15297acd-c406-46ae-ad52-27979f450abf
| Comm: tcp://127.0.0.1:43713 | Workers: 4 |
| Dashboard: /proxy/42597/status | Total threads: 16 |
| Started: Just now | Total memory: 44.92 GiB |
| Comm: tcp://127.0.0.1:37021 | Total threads: 4 |
| Dashboard: /proxy/44653/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:45285 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-w2qbeq7a | |
| Comm: tcp://127.0.0.1:43719 | Total threads: 4 |
| Dashboard: /proxy/34653/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:44203 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-5_jv0ml3 | |
| Comm: tcp://127.0.0.1:39265 | Total threads: 4 |
| Dashboard: /proxy/36927/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:39013 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-2rsbf9wf | |
| Comm: tcp://127.0.0.1:36493 | Total threads: 4 |
| Dashboard: /proxy/40415/status | Memory: 11.23 GiB |
| Nanny: tcp://127.0.0.1:44141 | |
| Local directory: /local/p93/rsp599/tmp/dask-worker-space/worker-5uyyp3ev | |
Wilma has kindly prepared the world ocean atlas data, currently 0.1$^\circ$ version available at /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/
! ls /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/
woa18_all_n01_10.nc woa18_all_o01_10.nc woa18_all_p01_10.nc woa18_all_n02_10.nc woa18_all_o02_10.nc woa18_all_p02_10.nc woa18_all_n03_10.nc woa18_all_o03_10.nc woa18_all_p03_10.nc woa18_all_n04_10.nc woa18_all_o04_10.nc woa18_all_p04_10.nc woa18_all_n05_10.nc woa18_all_o05_10.nc woa18_all_p05_10.nc woa18_all_n06_10.nc woa18_all_o06_10.nc woa18_all_p06_10.nc woa18_all_n07_10.nc woa18_all_o07_10.nc woa18_all_p07_10.nc woa18_all_n08_10.nc woa18_all_o08_10.nc woa18_all_p08_10.nc woa18_all_n09_10.nc woa18_all_o09_10.nc woa18_all_p09_10.nc woa18_all_n10_10.nc woa18_all_o10_10.nc woa18_all_p10_10.nc woa18_all_n11_10.nc woa18_all_o11_10.nc woa18_all_p11_10.nc woa18_all_n12_10.nc woa18_all_o12_10.nc woa18_all_p12_10.nc
!module load nco
!ncrcat /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_o*_10.nc /g/data/hh5/tmp/rsp599_tmp/nitrate.nc
ncrcat: ERROR no variables fit criteria for processing ncrcat: HINT Extraction list must contain a record variable to concatenate. A record variable is a variable defined with a record dimension. Often the record dimension, aka unlimited dimension, refers to time. To change an existing dimension from a fixed to a record dimensions see http://nco.sf.net/nco.html#mk_rec_dmn or to add a new record dimension to all variables see http://nco.sf.net/nco.html#ncecat_rnm
!ncdump -h /g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n01_10.nc
netcdf woa18_all_n01_10 {
dimensions:
xt_ocean = 3600 ;
yt_ocean = 2700 ;
st_ocean = 44 ;
time = 1 ;
variables:
double xt_ocean(xt_ocean) ;
xt_ocean:_FillValue = NaN ;
xt_ocean:long_name = "Nominal Longitude of T-cell center" ;
xt_ocean:units = "degree_east" ;
xt_ocean:modulo = 360. ;
xt_ocean:point_spacing = "even" ;
xt_ocean:axis = "X" ;
double yt_ocean(yt_ocean) ;
yt_ocean:_FillValue = NaN ;
yt_ocean:long_name = "tcell latitude" ;
yt_ocean:units = "degrees_N" ;
yt_ocean:cartesian_axis = "Y" ;
double st_ocean(st_ocean) ;
st_ocean:_FillValue = NaN ;
st_ocean:long_name = "zt" ;
st_ocean:units = "meters" ;
st_ocean:positive = "down" ;
st_ocean:point_spacing = "uneven" ;
st_ocean:axis = "Z" ;
double time(time) ;
time:_FillValue = NaN ;
time:long_name = "time" ;
time:cartesian_axis = "T" ;
time:calendar_type = "GREGORIAN" ;
time:units = "days since 0001-01-01" ;
time:calendar = "GREGORIAN" ;
double n_an(time, st_ocean, yt_ocean, xt_ocean) ;
n_an:_FillValue = NaN ;
n_an:long_name = "Objectively analyzed mean fields for moles_concentration_of_nitrate_in_sea_water at standard depth levels." ;
n_an:units = "micromoles_per_kilogram" ;
}
no3 = xr.open_mfdataset('/g/data/hh5/tmp/rsp599_tmp/nitrate_srf_mon_clim_2000_2019.nc')
no3
<xarray.Dataset>
Dimensions: (xt_ocean: 3600, yt_ocean: 2700, month: 12)
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
* month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
Data variables:
Nitrate (month, yt_ocean, xt_ocean) float32 dask.array<chunksize=(12, 2700, 3600), meta=np.ndarray>
Attributes:
long_name: nitrate
units: mmol/m^3
valid_range: [-1000000. 1000000.]
cell_methods: time: mean
time_avg_info: average_T1,average_T2...
QuantizeGranularBitRoundNumberOfSignificantDigits: 3
ncfiles: ['/g/data/cj50/access...
contact: Andrew Kiss
email: andrew.kiss@anu.edu.au
created: 2022-04-27
description: 0.1 degree ACCESS-OM2...
notes: Run configuration and...woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n01_10.nc')
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
<xarray.DataArray 'n_an' (st_ocean: 44, yt_ocean: 2700, xt_ocean: 3600)>
dask.array<getitem, shape=(44, 2700, 3600), dtype=float64, chunksize=(44, 2700, 3600), chunktype=numpy.ndarray>
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
* st_ocean (st_ocean) float64 0.5413 1.681 2.94 4.332 ... 630.7 695.4 766.1
Attributes:
long_name: Objectively analyzed mean fields for moles_concentration_of_n...
units: micromoles_per_kilogram# fig, axes = plt.subplots(ncols = 1, nrows=2, subplot_kw={'projection': ccrs.Robinson()}, figsize = (16, 6))
# cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
# 'fraction': 0.03,
# 'aspect': 15,
# 'shrink': 0.7}
# no3.Nitrate.isel(month=1).plot(ax=axes[0],
# x='xt_ocean', y='yt_ocean',
# transform=ccrs.PlateCarree(),
# vmin=0, vmax=30, extend='both',
# cmap=cm.cm.matter, cbar_kwargs=cbar_kwargs)
# woa_no3.n_an.isel(st_ocean=0).plot(ax=axes[1],
# x='xt_ocean', y='yt_ocean',
# transform=ccrs.PlateCarree(),
# vmin=0, vmax=30, extend='both',
# cmap=cm.cm.matter,
# cbar_kwargs=cbar_kwargs)
# axes[0].set_title('January Climatology ACCESS-OM')
# axes[1].set_title('January Climatology WOA18');
# import string
# # add features and label panels
# for n, ax in enumerate(axes):
# # features
# ax.coastlines(resolution='50m')
# ax.add_feature(land_50m)
# ax.gridlines(draw_labels=False)
# # subplot labels
# ax.text(0, 0.9, '('+string.ascii_lowercase[n]+')', transform=ax.transAxes, size=16)
no3diff = no3.Nitrate.isel(month=0) - woa_no3.n_an.isel(st_ocean=0)
no3diff.load()
<xarray.DataArray (yt_ocean: 2700, xt_ocean: 3600, time: 1)>
array([[[ nan],
[ nan],
[ nan],
...,
[ nan],
[ nan],
[ nan]],
[[ nan],
[ nan],
[ nan],
...,
[ nan],
[ nan],
[ nan]],
[[ nan],
[ nan],
[ nan],
...,
...
...,
[8.46696397],
[8.48541251],
[8.50233707]],
[[8.52319097],
[8.53394496],
[8.54389056],
...,
[8.47821153],
[8.49511876],
[8.51007693]],
[[8.51721955],
[8.52985752],
[8.54229768],
...,
[8.47746348],
[8.49057766],
[8.50500556]]])
Coordinates:
* xt_ocean (xt_ocean) float64 -179.9 -179.8 -179.7 ... 179.8 179.9 180.0
* yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 ... 89.89 89.94 89.98
month int64 1
st_ocean float64 0.5413
* time (time) object 0001-01-02 00:00:00fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Jan Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n02_10.nc')
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=1) - woa_no3.n_an.isel(st_ocean=0)
no3diff.load()
#
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Feb Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n03_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=2) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('March Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n04_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=3) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('April Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n05_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=4) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('May Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n06_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=5) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('June Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n07_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=6) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('July Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n08_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=7) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Aug Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n09_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=8) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Sep Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n10_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=9) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Oct Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n11_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=10) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Nov Clim (ACCES - WOA)');
woa_no3 = xr.open_mfdataset('/g/data/ik11/observations/woa18_bgc/woa18_bgc_10/woa18_all_n12_10.nc') # change filename
woa_no3 = woa_no3.assign_coords(xt_ocean=(((woa_no3.xt_ocean + 180) % 360) - 180)).sortby('xt_ocean')
woa_no3.n_an.squeeze(drop=True)
#
no3diff = no3.Nitrate.isel(month=11) - woa_no3.n_an.isel(st_ocean=0)### Change month name
no3diff.load()
# plotting
fig, axes = plt.subplots(ncols = 1, nrows=1, subplot_kw={'projection': ccrs.Robinson()}, figsize = (10, 6))
cbar_kwargs = {'label': '$NO_3 (mmolm^{-3}$)',
'fraction': 0.03,
'aspect': 15,
'shrink': 0.7}
axes.coastlines(resolution='50m')
axes.add_feature(land_50m)
axes.gridlines(draw_labels=False)
no3diff.plot(ax=axes,
x='xt_ocean', y='yt_ocean',
transform=ccrs.PlateCarree(),
robust=True, extend='both',
cmap=cm.cm.diff, cbar_kwargs=cbar_kwargs)
plt.title('Dec Clim (ACCES - WOA)');